Higher-order regression three-dimensional motion-compensation method for real-time optical coherence tomography volumetric imaging of the cornea

Abstract. Significance: Optical coherence tomography (OCT) allows high-resolution volumetric three-dimensional (3D) imaging of biological tissues in vivo. However, 3D-image acquisition can be time-consuming and often suffers from motion artifacts due to involuntary and physiological movements of the tissue, limiting the reproducibility of quantitative measurements. Aim: To achieve real-time 3D motion compensation for corneal tissue with high accuracy. Approach: We propose an OCT system for volumetric imaging of the cornea, capable of compensating both axial and lateral motion with micron-scale accuracy and millisecond-scale time consumption based on higher-order regression. Specifically, the system first scans three reference B-mode images along the C-axis before acquiring a standard C-mode image. The difference between the reference and volumetric images is compared using a surface-detection algorithm and higher-order polynomials to deduce 3D motion and remove motion-related artifacts. Results: System parameters are optimized, and performance is evaluated using both phantom and corneal (ex vivo) samples. An overall motion-artifact error of <4.61 microns and processing time of about 3.40 ms for each B-scan was achieved. Conclusions: Higher-order regression achieved effective and real-time compensation of 3D motion artifacts during corneal imaging. The approach can be expanded to 3D imaging of other ocular tissues. Implementing such motion-compensation strategies has the potential to improve the reliability of objective and quantitative information that can be extracted from volumetric OCT measurements.


Introduction
First introduced 30 years ago, optical coherence tomography (OCT) 1 has become a powerful medical imaging modality with micrometer resolution. It allows high-speed volumetric threedimensional (3D) imaging of biological tissues and has also been gaining a lot of interest as an intraoperative imaging modality to guide microsurgeries. [2][3][4] However, because a typical 3Dimage acquisition speed is in the order of a few hertz, it often suffers from motion artifacts due to involuntary and physiological movements of living tissue. 5 In particular, in ophthalmology, where OCT found its first uses and remains the primary field of application, 6 the target is constantly in motion. 7 This is because, even during so-called steady fixation on a target, there are involuntary fixational eye movements (e.g., tremors, drifts, and microsaccades). 8 Aside from such involuntary lateral ocular motion, another major source of involuntary eye movement is from vascular pulsation or respiration that manifests itself as motion in the axial direction 9 and may have a much more direct effect on the overall distortion of OCT images. 5 Any such involuntary eye movements may also be exaggerated in patients 10 and limit the reproducibility of quantitative measurements.
While a wealth of methods has been developed and explored to compensate for ocular motion during OCT imaging of the retina (both hardware-based 11,12 and software-based solutions [13][14][15][16][17][18] ), fewer works have focused on motion-compensated OCT-based measurements of the cornea. 5,19,20 Most of those are targeted toward corneal topography applications, 19,20 assume the surface to be perfectly spherical 20 or use a spherical model to estimate corneal curvature, 5 and are based on time-costly postprocessing algorithms. 5,18,19 In this paper, we propose and analyze a motion-compensated 3D-OCT approach geared toward real-time volumetric corneal imaging including in patients. Our approach is based on using a combination of intensity and topological analysis, where the tissue surface is identified and fitted to a polynomial curve and compared with three orthogonal references initially obtained to deduce the magnitude of motion-related image shift. The motion artifact along each (fast-axis) B-mode image is deduced by fitting these shift distances from the reference (slow-axis) images using high-order polynomials. This result is then used to compensate for the motion artifacts during the entire volumetric imaging procedure. To test and optimize system performance, we applied our algorithm to both corneal phantom and ex vivo bovine samples.

3D-OCT System
The 3D motion-compensated OCT system is based on an in-house built 100 kHz swept-source OCT (Axsun) system operating at the center wavelength of λ ¼ 1060 nm and a 110 nm tuning range. Figures 1(a) and 1(b) illustrate the OCT setup and proposed scanning protocol of the motion-compensated 3D-OCT system, respectively. The output of the swept laser is split into two arms by a 75:25 fiber coupler (TW1064R3A2A, Thorlabs). About 25% of the light is coupled to the sample arm and the other 75% is coupled to the reference arm. Polarization controllers (FPC020, Thorlabs) are used to adjust the polarization states of the two arms, enabling control of the interference between them. Achromatic collimators (AC300, Thorlabs) are utilized to collimate the light from the fiber. The collimated light is focused on the sample using an objective lens (LSM04, Thorlabs) with a numerical aperture (NA) of 0.05, corresponding to the lateral resolution of the OCT system of 0.61λ∕NA ¼ 12.9 μm. A vertical translation stage (VTS) (19089, Cenco-Lerner) allows simulation of axial sample motion (maximum around 1.5 cm) by precisely translating the stage using a micrometer. To simulate involuntary and random sample motion occurring in both lateral and axial dimensions, the sample was positioned in the palm of a hand. Galvanometer scanners (GVS002, Thorlabs), which contain two scanning mirrors are used to scan along the so-called "fast-axis" (X) and "slow-axis" (Y) to generate the volumetric data. The C-scan data composed of 1024 × 512 × 1024 voxels (X × Y × Z) corresponds to a 1 cm × 1.1 cm × 3.7 mm volume in corneal tissue. Note that the axial resolution, defined by the spectral bandwidth of the system, was experimentally confirmed to be 6 μm in air and 4.5 μm in corneal tissue. A dispersion compensation lens (LSM04DC, Thorlabs) is employed to compensate for the dispersion induced by the objective lens. Reference-and sample-arm signals are mixed by a 50:50 coupler (TW1064R3A2A, Thorlabs) and detected by a balanced detector. With future in vivo applications in mind, the power of the probe beam at the output (on the sample) was intentionally reduced to around 1 mW, to ensure light levels during imaging well below the maximum permissible exposure for ocular tissues. The signal-to-noiseratio at the imaging depth is >55 dB. The interferometric data is then processed by a data acquisition board (500 MSPS, 12-bit resolution) with a digitizer (embedded in the OCT engine from Axsun) before being transferred to a computer (Intel i7-9700F, Nvidia GeForce RTX 2060 SUPER, 16GB RAM) via a PCIe port by a frame grabber (PCIe-1433, National Instrument, Austin, Texas, United Statess) for further processing. Custom software written in C++/C# and a graphics processing unit (GPU) is used for high-speed data acquisition and processing, 3,21,22 resulting in a real-time B-scan rate of 100 Hz.
The proposed motion-compensated 3D-OCT system scanning protocol, as shown in Fig. 1(b), first scans three B-scans in the slow axis to obtain reference images, which are composed of 512 × 1024 pixels (Y × Z). This is immediately followed by a standard C-scan, scanning along X and Y, to collect volumetric images of the moving sample. The d is defined as the distance between the initial reference scans along the slow axis. L is the width of the B-scan image and α is defined as the ratio between d and L.

Workflow and Algorithm
After obtaining three B-scan reference images along the slow axis and as the C-scan data is being obtained, each B-scan image in the volumetric data is processed using a custom algorithm, first developed in MATLAB as a testbed prototype, and then programmed in C++ for real-time processing. The algorithm, as summarized in the flowchart depicted in Fig. 2, consists of two main threads, namely preprocessing and motion compensation, that are detailed in the following Secs. 2.2.1-2.2.3. Briefly, the preprocessing procedure is used to segment the surface boundary from the sample and fit to a higher-order polynomial curve. The motion-compensation procedure is based on this preprocessing result and uses it to deduce the sample motion by comparing it to the reference images. The effect of the motion is then subtracted from the image to display motion-free volumetric (C-scan) images. When multiple consecutive volumes are imaged, a new set of three references is acquired to achieve motion correction for each volume. Fig. 1 (a) A schematic representation of the OCT system setup. FC, fiber coupler; PC, optic fiber polarization controller; AC, achromatic collimator; GVS, galvanometer scanners that contain two scanning mirrors; OL, objective lens; VTS, vertical translation stage; DCL, dispersion compensation lens; BD+, BD-, balanced detector; DAQ, data acquisition board. (b) Schematic of the scanning protocol of the motion-compensated 3D-OCT system. d is the distance between adjacent reference images, L is the width of the B-scan image, and α is defined as the ratio between d and L.

Preprocessing algorithm
Before the thresholding can be applied, we first use a 10 × 10 Gaussian filter with a standard deviation of 4 for every B-scan image (I 0 ) to smooth the speckle pattern. Then, we remove the saturation artifacts introduced by light reflected from a highly specular surface that are over the dynamic range of the data-acquisition system. We simplified the saturation removal method based on the minimum variance mean-line subtraction described in Refs. 23 and 24, to reduce the computation cost. Specifically, I, the intensity value of the saturation-free image, is assigned as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 3 0 8 (1) where I g is the Gaussian-filtered image and A is the spatially averaging filter (size 40 × 1) with a value of 1/40. Next, to facilitate connected component analysis, we convert every grayscale B-scan into a binary image, based on a threshold value. More precisely, we applied the adaptive thresholding method 25 to calculate the thresholding value, and this ensures that the sample surface constitutes the largest area of bright components in the image after thresholding. Figure 3(a) shows the ex vivo bovine B-scan image. Figure 3(b) shows the corneal B-scan image after this thresholding was applied and clearly shows the upper and lower boundaries of the cornea sample.
The subsequent connected component analysis selects the regions with the largest connected elements of the maximal intensity in the image and erases the remaining bright and isolated pixels, retaining the basic shape of the sample boundaries, as illustrated in Fig. 3(c). To further improve this initial segmentation result, which is obtained from low-pass (Gaussian) filtered images, we also high-pass filtered the original B-scan images to obtain tissue boundaries. This result is used to further refine the segmented boundary result. Specifically, we first obtain high-pass filtered images (F) by applying a combination of Sobel filter, S ¼ ½1;2; 1; 0;0; 0; −1; −2; −1 T , and Prewitt filter, P ¼ ½1;1; 1; 0;0; 0; −1; −1; −1 T as follows: The preprocessing procedure is employed to segment the sample's surface boundary and fit to a higher-order polynomial. The motioncompensation procedure is based on this preprocessing result to deduce the sample motion by comparing it with the reference images, and then generates a motion-free C-scan by shifting every preprocessed B-scan accordingly.
Then, we search the high-pass filtered image, in the vicinity (AE15 pixels) of the initially segmented surface boundary [yellow line, Fig. 3(d)], for the pixel with the highest gradient (both horizontally and vertically) and use the result to replace the initial segmentation result. The resultant, corrected sample surface segmentation is fitted with a higher-order polynomial curve [red line, Fig. 3(d)]. The polynomial fitting algorithm uses the pixel location to form a Vandermonde matrix (V), resulting in a linear system VP ¼ Z, which is also Note that we define the points on the boundary with ðx i ; z i Þði ¼ 1; 2; : : : ; mÞ, where m is the number of points and n is the order of the polynomial. The coefficients of the polynomial are solved by P ¼ ðV T V −1 ÞV T Z, the matrix multiplication between the Moore-Penrose inverse of V and Z.

Lateral motion-compensation algorithm
The lateral motion-compensation algorithm is based on the preprocessed (segmented) result. We define the first B-scan cornea surface (z 1 ) as the target surface, then move every other B-scan surface (z t , t is the index of B-scan) horizontally to align it with the target surface. For each B-scan, the n'th order polynomial fitting curve of the corneal surface can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 1 8 5 The lateral movement of t'th B-scan could be figured out by minimizing cost function E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 4 0 in which, Δl t is the lateral movement, var represents variance, and m is the total number of A-scan (1024 in this case).

Axial motion-compensation algorithm
The axial motion-compensation algorithm uses the lateral motion-compensated result and is illustrated in Fig. 5. We first compare the B-scan surface detection result with the three reference images and obtain the differences in the axial position at the three intersecting planes of each B-scan. Then, the axial shift for every A-scan within the B-scan is estimated by fitting it with a polynomial and the shift is compensated accordingly in the final images. The motion-free volumetric image can be generated from sequentially-processed B-scans.
The three reference points [yellow dots in Fig. 5(a)] are defined as r i ðx i ; y i ; z i Þði ¼ 1; 2; 3Þ, which we write as r i ðx r i ; z r i Þði ¼ 1; 2; 3Þ, assuming all points in the same B-scan images share the  same y-axis coordinate. Then, we can deduce the axial motion at the three reference points, aðx r 1 Þ, aðx r 2 Þ, and aðx r 3 Þ as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 6 ; 7 1 1 The axial motion at every (other) points of zðxÞ, defined as mðxÞ is estimated by generating a second-order polynomial function based on the three reference points. Finally, we axially shift the location of every A-line in the B-scan to obtain the motion-compensated result [red line in Fig. 5(b)] fðxÞ as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 6 3 5 fðxÞ ¼ zðxÞ − aðxÞ: 3 Results We optimized system parameters and evaluated system performance mainly using phantom (thin plastic shell) and corneal (ex vivo bovine) samples. In addition, the order of polynomial for the surface fit was optimized using clinical images from a commercial ophthalmic OCT system (Avanti, Optovue) to assure it will work with various surface topographies. Tissue motion was simulated in two different ways. First, by positioning the samples on the vertical translation stage [VTS, Fig. 1(a)] and manually moving the stage up and down, axial motion was generated at frequencies ranging from about 0 to 5 Hz, with a typical amplitude around 1 mm. We know the frequency range of axial ocular motion could be from 0 to 125 Hz. 5 However, the amplitude of motion decreases rapidly with higher frequency and the axial resolution of our OCT system is 4.5 microns so that any smaller motions occurring at frequencies above about 5 Hz are not detectable with our system. Second, samples were positioned in the palm of the hand; such "free-hand" motion served to simulate the involuntary 3D motion of living tissue occurring with random frequencies and amplitudes.

Determining order of polynomial for corneal surface fit
With the proposed OCT system being geared toward ophthalmic applications in patients (e.g., to guide ocular surgeries and/or to extract quantitative information from volumetric measurements), our algorithm must work not only on normal corneas but also on pathological ones with abnormal surface topography. We hence used representative clinical OCT images (Avanti, Optovue), including from patients with bullous keratopathy and keratoconus to determine the optimal order of polynomial for the final surface fit during preprocessing of the data. Images were collected retrospectively; however, all subjects provided informed consent to have their images used in research, which was approved by the Institutional Review Board (Patient Protection Committee, Ile-de-France V) and adhered to the tenets of the Declaration of Helsinki. We applied different polynomial fits ranging from second-to sixth-order to the data (image) sets, during step (d) of the preprocessing algorithm [ Fig. 2(d)]. We used the fitting error, more precisely the root mean square error (RMSE) representing the difference between the surface fit and the ground truth, with the latter being estimated by manually segmenting the corneal surface, to quantify the results. The resulting fitting error for the three representative corneal images as a function of polynomial order is depicted in Figs. 6(d)-6(f). As expected, the fitting error decreases with increasing polynomial order for the healthy cornea with normal surface topography [ Fig. 6(d)]. However, this is not the case for the diseased corneas with abnormal surface pathology [see Figs. 6(e) and 6(f)], where one can note the error starting to increase again for the fifth-and sixth-order polynomials for the cornea with bullous keratopathy [ Fig. 6(e)] and keratoconus [ Fig. 6(f)] respectively. This may be due to the higherorder polynomials being very sensitive to the small, irregular surface variations of these diseased corneas. The fourth-order polynomial was hence chosen for the final surface fit in step (d) of the preprocessing procedure [see Fig. 2(d)].

Determining the number of reference planes
As detailed in Sec. 2.2.3, we deduce axial motion associated with every A-line within B-scan images with respect to three reference image planes. To arrive at three as the optimal number of reference planes for axial motion compensation, we tested our algorithm during C-scan imaging (consisting of 512 B-scans) of an ex vivo bovine cornea while inducing typical axial motion of around 1 Hz (such as induced from heartbeat) and changing the number of reference planes from 2 to 4. The reference planes are distributed based on the following rule, which minimizes the compensation error (as discussed in Sec. 3.1.3): suppose there are a total of t planes and L is the length of B-scan images, the i'th plane will be located at ð2i − 1ÞL∕2t. The time cost of the polynomial fitting process of the references [described by Eq. (3)] is measured and compared (using the built-in function in C++): Oðn 3 Þ, in which n is the dimension of the square matrix V in Eq. (3), and it is also the number of references in our case. In our test, first-, second-, and thirdorder polynomials were generated to estimate axial motion based on the 2, 3, and 4 reference points (planes), respectively. The RMSE between the motion-compensated and static (reference) C-scan image (calculated as the mean value of the absolute difference between the respective A-scan corneal surface positions) was used to quantify the results.
As summarized in Table 1, three references decrease the RMSE by around 4.21 μm with only an additional 2-ms time cost, compared with two references. Although four reference planes would even further reduce the RMSE by around 0.76 μm, it requires an additional 12 ms for processing and 0.76 μm is <1 pixel and the resolution of our imaging system. Hence, the three reference-plane approaches provide an optimum compromise between performance and time cost for our overall scanning speed and considered frequency range.

Determining the reference plane positions, α
Theory. The compensation result also depends on the position of the reference images and the distance between the adjacent reference images, d, and thus our selection of α, which is defined as d∕L, where L is the width of the raw B-scan image [see Fig. 1(b)]. Since we choose the middle reference to be fixed at the center of the x-direction, α can range from 0 (i.e., one reference plane at the center) to 0.5 (i.e., two references at the ends of the image). This is schematically shown in Figs. 7(a)-7(c), displaying the location of the three references (R1, R2, and R3; dotted lines) in a B-scan image, with α set as 0, 1/3, and 0.5, respectively. The quantitative analysis model is based on the hypothesis that the motion compensation accuracy based on references is a Gaussian-like function of the number of A-scans within the B-scan. The fitting error within every B-scan image is considered random, so the accuracy of  references used to compensate for the motion could be estimated as a Gaussian function (i.e., a widely used noise model to describe estimation-related error 26 ), and our experimental result agrees with this hypothesis. The accuracy function may be written as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 6 ; 6 8 7 in which, R is the accuracy, x is the number of A-scans, ranging from 0 to 1023 in our case, and i is the index for three references (i ¼ 1, 2, 3). We eliminated the parameters 1∕σ ffiffiffiffiffi 2π p to ensure the range of accuracy function is [0,1]. μ and σ are the mean value and the standard deviation of the Gaussian-like function, indicating the reference location and the width of the function shape, respectively. Therefore, μ of the three references can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 5 8 7 We assumed m ¼ 1023 to be the maximum number of A-scans. All three reference functions should have the same σ because all of them have the same compensation capability. Therefore, we have σ ¼ σ i ði ¼ 1; 2; 3Þ.
We define another accuracy function YðxÞ to determine the total accuracy of the motion compensation based on more than one reference. YðxÞ is a combination of several R i ðxÞ and should follow two rules: 1. When α → 0, YðxÞ → Rðx; n 2 ; σÞ; 2. When α → AE∞, YðxÞ → P i R i ðx; μ; σÞ.
These two rules ensure that YðxÞ also ranges from [0,1]. Thus, we construct YðxÞ as Yðx; μ i ; σÞ ¼ The number 1 in the denominator of distribution function ð1∕ðjx − μ i jÞ þ 1Þ is used to avoid zero. Figures 7(d)-7(f) show the accuracy functions R 1 ðxÞ, R 2 ðxÞ, R 3 ðxÞ, and YðxÞ with α set to 0, 1/3 and 0.5, respectively. Among the three examples, we obtain the best accuracy result with α ¼ 1∕3. The corresponding minimum accuracy of 0.33 represents the maximum error when the speed of the motion is approximately the same as the scanning rate of the reference images. However, since the reference images are obtained around 100 Hz, which is much faster than the anticipated sample motion, ranging up to 5 Hz, we expect our overall result to be highly accurate.
We also define an error function based on the accuracy function above to enable comparison with experimental data E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 2 4 2 Combining Eqs. (8) and (10) and assuming σ as a constant, we can write Eðμ i ; σÞ as a function of α E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 1 6 9 Eðα; σÞ ¼ X n x¼0 ½1 − Yðx; α; σÞ: (12) Experiment. To optimize α experimentally, we assessed the performance of the motion-compensation algorithm by performing several imaging trials while varying the distance between adjacent reference planes. The RMSE representing the difference between respective A-scan corneal surface positions of the motion-compensated and static (reference) C-scan image was used to quantify the results. Figure 8(a) shows the RMSE of the compensation result for around 1-Hz motion using two bovine corneal samples as a function of α. Two trials were performed for each sample. Figure 8(b) shows the RMSE result for different motion frequencies (0.4 to 5.2 Hz) using the corneal phantom as a function of α. Both experiments indicate that α around 0.33 achieves the best results (i.e., minimal RMSE) for all samples and frequencies used in the study. Note that, while the slower motion frequencies yielded smaller compensation errors (i.e., an RMSE of 2.34 for 1 Hz, and an RMSE of 1.26 for 0.4 Hz), even at higher frequencies (2.4 through 5.2 Hz), with α set to 0.33, our algorithm was capable of compensation with an error of <4.2 μm. Note that although the axial resolution of our system is 6 μm in air and 4.5 μm in corneal tissue, our ability to detect the peak of the point spread function is limited by the pixel size, which is 3.6 μm∕pixel (for corneal tissue). Hence, 4.2 μm corresponds to an overall compensation error of about 1 pixel.
Comparison. To further test the performance of our model, we calculated the error function with different σ [using Eq. (10)] and compared it with the experimental data, where the difference is the fitting error defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 3 6 7 WðσÞ ¼ with s and Expr being the index of α and the experimental data, respectively. As illustrated in Fig. 9(a), the fitting error (W) achieves a minimum at σ ¼ 123. Figure 9(b) depicts the result, at minimal σ ¼ 123, as a function of α (step size 10/1024), for the experimentally obtained motion compensation result and the numerical result for a motion frequency of 5.2 Hz; note that the data has been normalized. The average absolute error difference between  the theory and the experiment is 0.11. Both theory and experiment achieve a minimum error around α ¼ 0.338. Hence, we conclude that the mechanism of our motion compensation method is successfully described by our theoretical model.
The result above shows that the evenly-distributed reference planes can minimize the compensation error, which can be generalized to the distribution rule of reference planes we mentioned in Sec. 3.1.2.

Fourth-order polynomial for corneal surface fit
To validate the selection of the fourth-order as the optimal polynomial for the corneal surface fit, we assessed the performance of the algorithm experimentally. Specifically, two imaging trials using two ex vivo bovine corneal samples while inducing quasisinusoidal axial motion around 1 Hz were performed. The resulting images were fitted with the different order polynomials and compared. As previously done, the RMSE was used to quantify the results. The RMSE was deduced from measuring the difference between the reference and compensated C-scan image of the sample. The results are illustrated in Fig. 10 and show that fourth-order fitting achieved the best accuracy, i.e., the minimum RMSE. Considering this and our preliminary result using clinical images of human corneas with normal and abnormal surface topology, the fourth-order polynomial was confirmed as optimal surface fit in step (d) of the preprocessing procedure [ Fig. 2(d)].

Proof-of-concept demonstration of motion compensation
Using the optimized parameters, we performed a series of volumetric imaging trials using corneal phantom (thin plastic shell) and tissue samples (ex vivo bovine cornea), induced with artificially generated motions, as a proof-of-concept demonstration of the proposed motioncompensation method using our in house-built 3D-OCT system.
Axial motion. Figures 11(b), 11(e), 11(h), and 11(k) depict the volumetric (C-scan) images of the corneal phantom during induced axial motion [by means of the VTS, Fig. 1(a)] around 0.4, 2.4, 3.6, and 5.2 Hz, respectively, clearly showing the motion artifacts. Figures 11(a), 11(d), 11(g), and 11(j) are the corresponding static C-scan images, without any induced motion, that were used as references, and Figs. 11(c), 11(f), 11(i), and 11(l), depict the respective motion-compensated C-scan images using the proposed method. As expected, our algorithm successfully corrects the motion artifacts in these four motion frequency settings, with the smooth spherical surface of the phantom being visibly recovered throughout the entire volume of the sample. The RMSE, as quantitative metric of the results, is 1.3, 3.5, 3.7, and 4.2 μm for the axial motions of 0.4, 2.4, 3.6, and 5.2 Hz, respectively [see Fig. 8(b)].
An example of the axial motion-compensation result during volumetric imaging of an ex vivo bovine cornea is illustrated in Fig. 12. The average amplitude and frequency of the induced corneal motion were around 700 μm and 1 Hz, respectively. Figure 12(a) depicts the C-scan reference image of the cornea without any induced axial motion, whereas Figs. 12(b) and 12(c) are the C-scan images of the "moving" cornea before and after motion compensation. As expected, the motion-induced distortion of the corneal surface, clearly visible in Fig. 12(b), is successfully compensated using the proposed method [see Fig. 12(c)]. This example corresponds to trial 1 of sample 1 in Fig. 8(a), with an RMSE of 1.35 μm (for α ¼ 0.33). For direct visual comparison, the slow-axis cross-section showing the difference between reference and compensated result is also shown in the inset of Fig. 12(c).
3D motion. Figure 13 shows an example of the 3D motion-compensation result during C-scan (volume) "free-hand" imaging of a phantom [plastic shell; Figs   We also analyzed the lateral and axial motion that was induced during free-hand imaging of the corneal samples. The amplitudes of axial and lateral motion occurring during free-hand imaging of the ex vivo bovine cornea from Fig. 13 are depicted in Figs. 14(a) and 14(b), respectively; the corresponding frequencies are shown in Fig. 14(c). The mean value of axial and lateral motion between adjacent B-scan images is around 14.18 and 165 μm, respectively. Also note that most energy of motion, in either direction, occurs in the range from 0 to 5 Hz and that motion above 15 Hz is negligible.

Algorithm Efficiency
Our algorithm was first developed in MATLAB as a testbed prototype, then implemented in C++ for real-time processing. To demonstrate its efficiency and real-time capabilities, we evaluated the time consumption of the algorithm during volumetric (C-scan) imaging. The tested C-scan image contains 512 B-scans and the size of each B-scan is 1024 × 1024. The mean value of time consumption for each B-scan is 7.6 ms, including both preprocessing and motion compensation procedures, which is faster than our 100-Hz B-scan rate.
Detailed time consumptions for the diverse processing steps are shown in Fig. 15. Approximately 36% of the 7.6 ms are used for thresholding, 28% are used for various filter operations based on convolution, 26% are used for connect component analysis, 9% are used for matrix-related calculations, and the remaining 1% are used for other procedures.  The preprocessing procedure is much more time-consuming than the actual motion compensation because the latter is based on a single variable polynomial, which only requires few computational resources.
To further improve algorithm efficiency, the GPU driven by compute unified device architecture is assigned to the thresholding task, matrix calculation, and filter operation, which are mainly responsible for the processing time. This parallel processing reduces the processing time considerably, resulting in a mean value of overall time consumption for each B-scan image preprocessing and motion compensation of 3.40 ms, including data I/O between CPU and GPU.
A single B-scan image is processed simultaneously during the acquisition of the subsequent B-scan frame. With the processing time (3.40 ms) being much smaller than the actual acquisition time (10 ms, set by the OCT engine), every motion-compensated B-mode frame is processed and displayed in real-time, at the same speed as the image acquisition.

Discussion
We implemented an OCT system for volumetric imaging of the cornea, capable of compensating both axial and lateral motion using a combination of intensity and topological analysis, with the corneal surface being identified and fitted to a polynomial curve and compared with three orthogonal references (slow-axis) B-scans initially obtained to deduce the magnitude of motionrelated shifts to correct the fast-axis B-cans throughout the OCT volume.
The use of orthogonal B-scans to compensate for motion in OCT volumes is not new, especially for the retina. For example, Kraus et al. 13 retrospectively gathered additional B-scans along the slow axis of the acquired OCT volume to align the fast-axis (orthogonal) B-scans, on a per A-scan basis. Other related retinal methods include those by Lezama et al., 27 who also used a segmentation approach, and by Kim et al. 18 who performed simultaneous orthogonal scanning by means of two perpendicularly polarized beams. Our method could also be applicable for retinal imaging, including OCT angiography (OCT-A), upon redesigning the scanning and microscope objective and imaging system of the OCT apparatus. Note that due to the refractive power of the cornea and crystalline lens, the OCT image of the retina is distorted and does not represent the true curvature of the retina so that additional image processing is required to obtain the true reference image. However, for many applications (including robotic tool guidance that we are pursuing 2,28-30 ), an optically distorted image is sufficient and only the relative position and the curvature are needed to correctly assess the tool to tissue distance. Therefore, even for retinal imaging, one could still use the same processing methods to obtain motion-free retinal images.
An advantage of our approach is that every target surface is represented as a polynomial curve, so the difference between different B-scans can be calculated with low computational cost by the comparison of polynomial parameters rather than complicated and time-consuming cross-correlation, 18,27 or local pairwise phase correlation. 31 Although the simulated motion in this preliminary (ex-vivo) proof-of-concept study does not reflect the complex movements of the living eye, most involuntary and physiological ocular (lateral and axial) motion is expected to fall within the range from 0 to 5 Hz (including micro saccades 32 ), which is much smaller than our B-scan rate of 100 Hz. Hence, sample movement within each B-scan image could be approximated by a second-order polynomial curve, only requiring three references and reducing scanning time cost.
The use of higher-order polynomial, rather than spherical approximations, 5 for corneal surface fits, also ensures the recognition and differentiation of involuntary tissue movement from abnormal surface topology, which may be particularly important for imaging in patients with corneal disease.
Similarly, with one of the goals of developing the algorithm has been to use the motioncompensated imagine system for ophthalmic applications in patients (e.g., to guide the robotics eye surgery and/or to extract quantitative information from volumetric measurements), a raster scan was favored over a radial-scan approach for its ease of use (e.g., to locate and track the movement of the surgical tool, such as needles).
Limitations of the current approach include the assumption of no motion occurring during the acquisition of the three reference (slow-axis) B-scans. In fact, the observed increase in RMSE with higher-frequency motion [ Fig. 8(b)] could be the result of higher-velocity movements affecting the references and may become enhanced during in-vivo application of the approach, in particular during imaging of patients whose involuntary eye motion may be exaggerated. 10 To reduce potential motion artifacts during reference acquisition, we plan to decrease the number of A-scans during the acquisition of the references (slow-axis B-scan images) during future in vivo studies. For example, the current size of a reference B-scan is 512 × 1024, with an acquisition time of 5 ms, which is shorter than the typical (fast-axis) B-scan image acquisition time of 10 ms. However, further reducing the reference B-scans to only 100 A-scans would still be sufficient to perform the motion compensation using our algorithm and reducing the acquisition time for each reference to around 1 ms so that motion during reference-plane acquisition could indeed be considered negligible.
The algorithm could also be accelerated by choosing a smaller region of interest (ROI). Our C-scan data composed of 1024 × 512 × 1024 voxels (X × Y × Z) corresponds to a 1 cm × 1.1 cm × 3.7-mm corneal tissue volume. If we reduced the ROI to around 0.5 cm × 0.65 cm × 3.7 mm with the same image resolution, e.g., we would only spend one-fourth time for the image acquisition and processing and 1/2 time for reference acquisition. This would also reduce the effect of motion, in particular higher-frequency motion that might be pronounced during patient (in-vivo) imaging, and lead to higher compensation accuracy.
Future work will include such optimizations to address these limitations, enhance the approach and optical system to be more clinically oriented (including accounting for more complex motion of the living eye), and test system performance during in vivo measurements.

Conclusion
We described a 3D motion-compensated OCT system based on higher-order regression toward real-time volumetric imaging of the cornea. Fourth-order polynomial surface fitting in combination with the acquisition of three orthogonal (slow-axis) reference planes ensures differentiation of abnormal surface topology from involuntary tissue motion. Volumetric (C-scan) imaging of corneal phantom and bovine samples during simulated (axial and free-hand) motion demonstrated effective compensation of 3D motion artifacts with an overall error of <4.61 μm. The processing time for each (fast-axis) B-scan is 3.40 ms using a GPU-based parallel computation technique. Our approach can be expanded to 3D-imaging of other ocular tissues. Implementing such motion-compensation strategies into clinical OCT systems has the potential to improve the reliability of objective and quantitative information that can be extracted from volumetric OCT measurements.

Disclosures
The authors declare that there are no conflicts of interest related to this article.